% Script for creating speed profile plots

%% Options

clear
clc

file = 'H:\My Drive\Boats\ReplicationCode\data\final\all_points' ;
out_dir = "H:\My Drive\Boats\ReplicationCode\results\voyage_profiles" ;

%% Figure Options %%%%%%%%%%%%%%%%
font_size = 12 ;
figure_line_width = 2.5 ;
figure_axis_width = 1.5 ;
legLineWidth = 0.5 ;

%%% Change Default Color Map
%cmap = linspecer(6,'qualitative');  % colormap('colorcube')

%%% colors consistent with stata scheme
cmap_stata = [ 85 117 47
    227 126 0
    26 71 111
    144 53 59
    110 142 132
    193 5 52 ] / 255 ; 


%%% Setup for output
%%% set defaults 
set(0, 'defaultAxesColorOrder',cmap_stata)
set(0, 'DefaultTextFontSize' , font_size )
set(0, 'DefaultAxesFontSize' , font_size )
set(0, 'defaultlinelinewidth', figure_line_width)
set(0, 'defaultaxeslinewidth', figure_axis_width)
set(0, 'DefaultAxesTickDir', 'out')
set(0, 'DefaultAxesTickLength', [0.01,0.025] )
set(0, 'DefaultFigureWindowStyle', 'docked') %%% making added figures be docked

set(groot, 'DefaultTextInterpreter', 'latex')
set(groot, 'DefaultLegendInterpreter', 'latex')

% size of figures
fig_height = 7 ;
fig_width = fig_height*(4/3) ; % making 4:3 aspect ratio

markers = {'o','d','^','s'} ;

%% Set up
if ~exist(out_dir, 'dir')
   mkdir(out_dir)
end

%% Loading and cleaning

disp('Loading')
%Tfull = readtable(file,'DatetimeType','text') ;
Tfull = table() ;
for y=2009:2016 
    disp(y)
    Tfull = [ Tfull ; readtable( [ file,'_',mat2str(y),'.csv'],'DatetimeType','text')  ] ;
end


disp('Processing')
% Need to update OID (could be duplicate OID over years)
Tfull.oid_old = Tfull.oid ;
Tfull.oid = findgroups(Tfull.oid,Tfull.aisyear) ;

% convert to times
Tfull.time1 = datetime( Tfull.time1 , 'InputFormat', 'yyyy-MM-dd HH:mm:ss.S') ;
Tfull.time2 = datetime( Tfull.time2 , 'InputFormat', 'yyyy-MM-dd HH:mm:ss.S') ;
Tfull.date_ds_apeep = datetime( Tfull.date_ds_apeep , 'InputFormat', 'yyyy-MM-dd HH:mm:ss') ;

%{
* generating year, month, and month*year indicators
%}
Tfull.WC1 = Tfull.port1<=12 ;
Tfull.WC2 = Tfull.port2<=12 ;

Tfull.time = Tfull.time1 ;
Tfull(Tfull.WC1==0 & Tfull.WC2==0,'time') = Tfull(Tfull.WC1==0 & Tfull.WC2==0,'time2') ;

ecaCA1 = datetime(2009,7,1) ; % established
ecaCA2 = datetime(2011,12,1) ; % boundary change
ecaNA1 = datetime(2012,8,1) ; % established
ecaNA3 = datetime(2015,1,1) ; % goes to 0.1%

eca_names = {'Pre ECA','CA 2009','CA 2011','CA+NA (1)','CA+NA (0.1)'};

Tfull.ecaCA1 = Tfull.time >= ecaCA1 ;
Tfull.ecaCA2 = Tfull.time >= ecaCA2 ;
Tfull.ecaNA1 = Tfull.time > ecaNA1 ;
Tfull.ecaNA3 = Tfull.time > ecaNA3 ;
Tfull.eca_ind = Tfull.ecaCA1 + Tfull.ecaCA2 + Tfull.ecaNA1 + Tfull.ecaNA3 ; 

Tfull.eca_name = eca_names(Tfull.eca_ind+1)';

% flagging steam turbines
Tfull.steam = Tfull.powertype == "Steam Turbine" ;

% finding exempt vessels
Tfull.exempt = Tfull.steam | ( Tfull.enginedispl_cul<30 & Tfull.loa<121 & Tfull.gt<10000) ;
Tfull.treat = ~Tfull.exempt ;

% update nan to 0 here, because these are nan for first segment
Tfull(isnan(Tfull.in_eca09),'in_eca09')= {0} ;
Tfull(isnan(Tfull.in_eca11),'in_eca11')= {0} ;

Tfull.in_eca = Tfull.in_eca09 .* (Tfull.eca_ind<=1) + Tfull.in_eca11 .* (Tfull.eca_ind>1) ;

% number of unique OID by route, and ECA treatment
Troute_count = varfun(@(x) numel(unique(x)),Tfull,'GroupingVariables',{'port1','port2','eca_name','group_agg'},'InputVariables','oid') ;

% 
%Tfull.m = num2str( month(Tfull.TIME1) ) ;
Tfull.my = findgroups( datenum( year(Tfull.time) , month(Tfull.time) , 1 ) ) ;
Tfull.t = findgroups( datenum(year(Tfull.time) , month(Tfull.time) , day(Tfull.time) ) ) ;

Tfull.t_eca1 = datenum(year(Tfull.time) , month(Tfull.time) , day(Tfull.time) ) - datenum(ecaCA1) ;
Tfull.t_eca2 = datenum(year(Tfull.time) , month(Tfull.time) , day(Tfull.time) ) - datenum(ecaCA2) ;
Tfull.t_eca3 = datenum(year(Tfull.time) , month(Tfull.time) , day(Tfull.time) ) - datenum(ecaNA1) ;
Tfull.t_eca5 = datenum(year(Tfull.time) , month(Tfull.time) , day(Tfull.time) ) - datenum(ecaNA3) ;

Tfull.SBind = cellfun( @(x) ~isempty(x), strfind(Tfull.port_id,'51') )  ;

Tfull.sample_wexempt = Tfull.kmh<100 & Tfull.km_onland<5 & ( Tfull.ensenada==0 & Tfull.rosarito==0 & Tfull.elsegundo==0 ) ;
Tfull.sample = Tfull.sample_wexempt & Tfull.treat==1 ;


% getting cumulative distance in reverse direction
Tfull = sortrows(Tfull,{'oid','date_ds_apeep'}) ;
a = varfun(@(x) cumsum(x,1,'reverse') ,  Tfull , 'GroupingVariables' , 'oid' , 'InputVariables' , 'd_km' ) ;
Tfull.km_cum_rev = a.Fun_d_km ;

Tfull.km_cum_adj = Tfull.km_cum .* ( Tfull.dirnw ) + Tfull.km_cum_rev .* ( 1 - Tfull.dirnw  ) ;
Tfull.km_cum_adj_naeca = Tfull.km_cum_adj .* ( Tfull.in_naeca>0) ;
a = varfun(@max,Tfull,'GroupingVariables','oid','InputVariables','km_cum_adj_naeca') ; 
Tfull = innerjoin(Tfull,a,'Keys','oid','RightVariables','max_km_cum_adj_naeca') ;
clear a 
Tfull.km_to_naeca = Tfull.km_cum_adj - Tfull.max_km_cum_adj_naeca  ;

Tfull(:,'dum_ones')= {1} ;

%% Average speed and share within ECA
routes = { 'LALB_SFBay' , [3,4] , [1]  , 400 
           'SFBay_Sea'  , 1 , 7 , 400 
           'SFBay_West'  , [1] , [103,104,105,106] , 200 
           'LALB_West'  , [3,4] ,[103,104,105,106] , 400
          } ;

vesTypes = {
    'Container',{'FCC'} , [18,42] ...
    %'Bulker',{'Bulker'} , [10,30] ...
    %'OtherCargo',{'Other Cargo'} , [10,30] ...
   % 'Tanker',{'Tankers'} , [10,30]
    } ;

pols_plot = {'PreECA','CA2009','CA2011'} ; 
        
for v = 1:size(vesTypes,1)        
    for r = 1:size(routes,1) 
        port1 = routes{r,2} ;
        port2 = routes{r,3} ;
        T = Tfull( ismember(Tfull.port1,port1) & ismember(Tfull.port2,port2) & ismember(Tfull.group_agg,vesTypes{v,2}) & Tfull.sample==1 , :) ;

        ylims = vesTypes{v,3} ;
        ytick = ylims(1):4:ylims(2) ;
        bw=[0,routes{r,4}];
        bin = 20 ;
        bins = bw(1):bin:bw(2) ;
        T.km_bins = discretize(T.km_cum,bins,bins(2:end)) ;

        Tbins = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name'},'InputVariables',{'kmh'}) ;
        Tbins.eca_name = cat(2 , char( Tbins.eca_name )   ) ;
        % this breaks apart by SB IND
        %Tbins = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name','SBind'},'InputVariables',{'kmh'}) ;
        %Tbins.eca_name = cat(2 , char( Tbins.eca_name ) ,   num2str( Tbins.SBind )   ) ;
        %Tbins.SBind = [] ;
        Tbins.GroupCount = [] ;
        Tbins = unstack(Tbins,{'mean_kmh'},'eca_name') ;
        
        Tbins_ineca = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name'},'InputVariables',{'in_eca'}) ;
        Tbins_ineca.eca_name = cat(2 , char( Tbins_ineca.eca_name )    ) ;        
%         Tbins_ineca = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name','SBind'},'InputVariables',{'in_eca'}) ;
%         Tbins_ineca.eca_name = cat(2 , char( Tbins_ineca.eca_name ) ,   num2str( Tbins_ineca.SBind )   ) ;
%         Tbins_ineca.SBind = [] ;
        Tbins_ineca.GroupCount = [] ;
        Tbins_ineca = unstack(Tbins_ineca,{'mean_in_eca'},'eca_name') ;
        
        Tbins(:,pols_plot(~ismember(pols_plot,Tbins.Properties.VariableNames))) = {NaN};
        Tbins_ineca(:,pols_plot(~ismember(pols_plot,Tbins_ineca.Properties.VariableNames))) = {NaN};
        
        track_cnt = numel(unique(T.oid)) ;
        
        fig = figure() ;    
        %[ax,h1,h2] = plotyy(Tbins.km_bins,Tbins{:,{'PreECA1','CA20090','CA20091'}},Tbins_ineca.km_bins,Tbins_ineca{:,{'PreECA1','CA20090','CA20091'}},'plot') ;
        [ax,h1,h2] = plotyy(Tbins.km_bins,Tbins{:,{'PreECA','CA2009','CA2011'}},Tbins_ineca.km_bins,Tbins_ineca{:,{'PreECA','CA2009','CA2011'}},'plot') ;
        set(h2, {'Color'}, get(h1,'Color'));
        set(h2, 'LineStyle' , '--' ) ;
        
        for i=1:numel(h1)
            set(h1(i),'Marker' , markers{i} );
            set(h2(i),'Marker' , markers{i} );
        end

        lw = get(h2,'LineWidth');
        set(h2, 'LineWidth' , lw{1}/2 ) ;
        ylabel(ax(1),'Speed (km/hr) - Solid')
        ylabel(ax(2),'Share in ECA - Dashed')
        xlabel('Distance (km)')
        ylim(ax(1),ylims)
        yticks(ax(1),ytick)
        yticks(ax(2),0:0.25:1)
        legend(h1,eca_names,'Orientation','Horizontal','Location','southoutside')
        set(gca,'box','off')

        [fileName] = printFig(fig, ['speed_ineca_',vesTypes{v,1},'_',routes{r,1}] , out_dir , fig_width , fig_height, 0) ;

    end
end

%% Speed profiles of treated vs untreated vessels 

routes = { 'SFBay_Sea'  , 7 , 1 , 300  , 1
           'SFBay_West'  , [1] , [105] , 200 , 0
           'LALB_West'  , [3,4] , [103] , 300 , 0
          } ;

vesTypes = {
    'Container',{'FCC'}  , [18,42]  ... [18,42]
    } ;

pols_plot = {'PreECA','CA2009','CA2011'} ; 
        
for v = 1:size(vesTypes,1)  
    for treat = [0,1] % looping through exempt and treated
        for r = 1:size(routes,1) 
            port1 = routes{r,2} ;
            port2 = routes{r,3} ;
            T = Tfull( ismember(Tfull.port1,port1) & ismember(Tfull.port2,port2) & ismember(Tfull.group_agg,vesTypes{v,2}) & Tfull.sample_wexempt==1 & Tfull.treat==treat , :) ;

            ylims = vesTypes{v,3} ;
            ytick = ylims(1):4:ylims(2) ;
            bw=[0,routes{r,4}];
            bin = 20 ;
            bins = bw(1):bin:bw(2) ;
            if routes{r,5}==0
                T.km_bins = discretize(T.km_cum,bins,bins(2:end)) ;
            else
                T.km_bins = discretize(T.km_cum_adj,bins,bins(2:end)) ;
            end

            var = 'kmh' ; % 'distuscoast_km'
            Tbins = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name'},'InputVariables',{var}) ; % kmh
            Tbins.eca_name = cat(2 , char( Tbins.eca_name )   ) ;
            % this breaks apart by SB IND
            %Tbins = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name','SBind'},'InputVariables',{'kmh'}) ;
            %Tbins.eca_name = cat(2 , char( Tbins.eca_name ) ,   num2str( Tbins.SBind )   ) ;
            %Tbins.SBind = [] ;
            Tbins.GroupCount = [] ;
            Tbins = unstack(Tbins,{['mean_',var]},'eca_name') ;

            Tbins_ineca = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name'},'InputVariables',{'in_eca'}) ;
            Tbins_ineca.eca_name = cat(2 , char( Tbins_ineca.eca_name )    ) ;        
    %         Tbins_ineca = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name','SBind'},'InputVariables',{'in_eca'}) ;
    %         Tbins_ineca.eca_name = cat(2 , char( Tbins_ineca.eca_name ) ,   num2str( Tbins_ineca.SBind )   ) ;
    %         Tbins_ineca.SBind = [] ;
            Tbins_ineca.GroupCount = [] ;
            Tbins_ineca = unstack(Tbins_ineca,{'mean_in_eca'},'eca_name') ;

            Tbins(:,pols_plot(~ismember(pols_plot,Tbins.Properties.VariableNames))) = {NaN};
            Tbins_ineca(:,pols_plot(~ismember(pols_plot,Tbins_ineca.Properties.VariableNames))) = {NaN};

            track_cnt = numel(unique(T.oid)) ;

            fig = figure() ;    
            %[ax,h1,h2] = plotyy(Tbins.km_bins,Tbins{:,{'PreECA1','CA20090','CA20091'}},Tbins_ineca.km_bins,Tbins_ineca{:,{'PreECA1','CA20090','CA20091'}},'plot') ;
            [ax,h1,h2] = plotyy(Tbins.km_bins,Tbins{:,{'PreECA','CA2009','CA2011'}},Tbins_ineca.km_bins,Tbins_ineca{:,{'PreECA','CA2009','CA2011'}},'plot') ;
            set(h2, {'Color'}, get(h1,'Color'));
            set(h2, 'LineStyle' , '--' ) ;

            for i=1:numel(h1)
                set(h1(i),'Marker' , markers{i} );
                set(h2(i),'Marker' , markers{i} );
            end

            lw = get(h2,'LineWidth');
            set(h2, 'LineWidth' , lw{1}/2 ) ;
            ylabel(ax(1),'Speed (km/hr) - Solid')
            ylabel(ax(2),'Share in ECA - Dashed')
            xlabel('Distance (km)')
            ylim(ax(1),ylims)
            yticks(ax(1),ytick)
            yticks(ax(2),0:0.25:1)
            legend(h1,eca_names,'Orientation','Horizontal','Location','southoutside')
            set(gca,'box','off')


            [fileName] = printFig(fig, ['speed_ineca_',vesTypes{v,1},'_',routes{r,1},'_treat',num2str(treat)] , out_dir , fig_width , fig_height, 0) ;

        end
    end
end


%% Marginal Damage Profiles

% name, port ids orig, port ids dest  
routes = {  'LALB_SFBay' , [3,4] , [1]  
          } ;

pol = 'pm_inm' ;
label = 'Marginal Damage of PM (\$1000/ton)' ;
scale = 1000 ;
yrange = [0,400] ;

Tfull.D_plot = Tfull.(['d_',pol]) / scale ;

for r = 1:size(routes,1) 
    port1 = routes{r,2} ;
    port2 = routes{r,3} ;
    T = Tfull( ismember(Tfull.port1,port1) & ismember(Tfull.port2,port2) & Tfull.steam==0 & ismember(Tfull.group_agg,'FCC') & Tfull.kmh<100 , :) ;
    
    bw=[0,250];
    bin = 25 ;
    bins = bw(1):bin:bw(2) ;
    T.km_bins = discretize(T.km_cum,bins,bins(2:end)) ;
    Tbins = varfun(@mean,T,'GroupingVariables',{'km_bins','eca_name'},'InputVariables',{'D_plot'}) ;
    Tbins.GroupCount = [] ;
    Tbins = unstack(Tbins,{'mean_D_plot'},'eca_name') ;

    fig=figure()
    h=plot(Tbins.km_bins,Tbins{:,{'PreECA','CA2009','CA2011'}}) ; % 
    for i=1:numel(h)
        set(h(i),'Marker' , markers{i} );
    end
    
    ylabel(label)
    ylim(yrange)
    xlabel('Distance (km)')
    legend(h,eca_names,'Orientation','Horizontal','Location','southoutside')
    set(gca,'box','off')
    [fileName] = printFig(fig, ['D_PM_prof_',routes{r,1}] , out_dir , fig_width , fig_height, 0) ;

end
